Mapping genotype to phenotype information across magnifications¶
In [43]:
import tifffile as tf
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
from scipy import spatial
import pandas as pd
import seaborn as sns
import PIL
from math import sqrt
import matplotlib.patheffects as PathEffects
import random
from tqdm import tqdm
from scipy.optimize import minimize
from PIL import Image
from os import listdir
from os.path import isfile, join
import warnings
In [44]:
from helper_functions import *
In [46]:
#this alters the amount of width taken up but notebook chunks
from IPython.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))
Set well that will be mapped between magnifications¶
In [47]:
well_ID = 12
Retrieve coordinates associated with tile centers at each magnification¶
In [48]:
lomag_tiles_dir = '/path/to/your/raw/low_magnification/tifs'
himag_tiles_dir = '/path/to/your/raw/high_magnification/tifs'
In [49]:
#replace these with your paths at runtime
lomag_tiles_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/input_files/raw/genotyping_tifs/cycle_0/'
himag_tiles_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/input_files/raw/phenotyping_tifs/Mito_wells/'
In [50]:
lomag_df = RetrieveCoordinates_in_df(lomag_tiles_dir, well_ID)
himag_df = RetrieveCoordinates_in_df(himag_tiles_dir, well_ID)
In [51]:
lomag_df.head()
Out[51]:
Use tile center coordinates from each magnification to identify pairs of tiles that correspond to approximately the same location in the well at each magnification¶
In [52]:
#calculate radius and centroid of higher magnification tiles
_well_radius, _well_centroid = ApproximateRadiusAndCentroid(himag_df)
In [53]:
well_fiducial_tiles = GenerateFiducials_df(lomag_df, himag_df, _well_radius, _well_centroid)
well_fiducial_tiles.head()
Out[53]:
Use tile IDs to look up images at each magnification and search for the same cell at the second magnification¶
In [55]:
#both dataframes are generated from OPS outputs from processing pipelines documented in this GitHub repo: feldman4/OpticalPooledScreens
#replace these with your output at runtime
df_lomag = pd.read_csv('OPS_genotyping_output.csv', index_col=0) #readcalling output from OPS pipeline
df_himag = pd.read_csv('OPS_phenotyping_output.csv') #output from OPS pipeline
In [56]:
df_lomag.sample(3)
Out[56]:
In [57]:
df_himag.sample(3)
Out[57]:
In [58]:
#outputs from standard OPS processing of raw images
lomag_processed_images_dir = '/path/to/your/OPS_processed/low_magnification/tifs'
himag_processed_images_dir = '/path/to/your/OPS_processed/high_magnification/tifs'
In [59]:
#replace these with your paths at runtime
lomag_processed_images_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/Readcalling/Mito_process_05072023_thresh90_v1/'
himag_processed_images_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/Phenotyping/Mito_process_05072023_pheno_v1/'
In [60]:
#instantiate empty lists to fill with local coordinate information
cellpairs_list = []
tilepairs_list = []
In [63]:
#for each pair of tiles in the well_fiducial_tiles dataframe, print out images and identify the number associated with the same cell at each magnification
PrintImagesWithAllNucCentroids(tile_lomag=157, tile_himag=2385, lomag_image_dir=lomag_processed_images_dir, himag_image_dir=himag_processed_images_dir,
cells_lomag_df=df_lomag, cells_himag_df=df_himag,
well=well_ID, suppress_labels=False, sample_lomag=0.8, _vmax=2500)
Out[63]:
In [64]:
#after investigation of above images, supply manually identified fiducial cell identities
#cell_IDs
cellpair_lomag_himag = [174,12] #of form [cell_ID_in_lomag_image, cell_ID_in_himag_image]
cellpairs_list.append(cellpair_lomag_himag)
#tile_IDs - these will typically match the tileIDs provided in the dataframe, but on occaision, you may choose to search neighboring tiles if a shared cell can't be identified
tilepair_lomag_himag = [157,2385]
tilepairs_list.append(tilepair_lomag_himag)
In [65]:
#repeat this for every line in the well_fiducial_tiles dataframe
#after repeating this for 8-12 fiducials total:
cellpairs_list = [[1683, 36], [160, 8], [461, 19], [863, 22], [2587, 97], [1911, 92], [17, 16], [129, 20], [2063, 62], [45, 8], [675, 22], [655, 12]]
tilepairs_list = [[8, 64], [157, 2385], [89, 1227], [76, 1171], [32, 298], [23, 334], [133, 2115], [142, 2148], [40, 499], [125, 1947], [86, 1238], [79, 1182]]
In [66]:
fiducial_cells = pd.DataFrame({
'cell_lomag': [pair[0] for pair in cellpairs_list],
'cell_himag': [pair[1] for pair in cellpairs_list],
'tile_lomag': [pair[0] for pair in tilepairs_list],
'tile_himag': [pair[1] for pair in tilepairs_list],
'well': well_ID})
In [67]:
#look up tilewise coordinates for each of these cells to later place on global coordinate scale
points_hi, points_lo = GetTilewiseLocations(df=fiducial_cells, himag_df=df_himag, lomag_df=df_lomag, well=well_ID)
points_hi, points_lo
Out[67]:
Use the coordinates of the fiducial cells to generate a translation, rotation, scaling matrix between the magnifications¶
In [68]:
# Get image height and width
example_tif_lomag = os.path.join(lomag_tiles_dir, random.choice([f for f in os.listdir(lomag_tiles_dir) if f.endswith('.tif')]))
example_tif_himag = os.path.join(himag_tiles_dir, random.choice([f for f in os.listdir(himag_tiles_dir) if f.endswith('.tif')]))
h_lomag, w_lomag = Load_TIF(example_tif_lomag).shape
h_himag, w_himag = Load_TIF(example_tif_himag).shape
In [69]:
# Generate a starting grid to represent how tiles relate to each other in 2d space.
# This will depend on the manner in which the stage moves as images are collected, and will be investigated later and corrected if necessary
untested_grid_lomag, lomag_coord_df = GenerateTileGrid(path_name=lomag_tiles_dir, well=12)
untested_grid_himag, himag_coord_df = GenerateTileGrid(path_name=himag_tiles_dir, well=12)
In [70]:
print(untested_grid_lomag)
In [71]:
#let's test arrangment by picking a set of four tiles that are adjacent in the grid, and determining which transformation (if any) is required to align them
input_dict = {'top_left':81,'top_right':82,'bottom_left':98,'bottom_right':97}
TestTileGrid(input_dict, lomag_tiles_dir, lomag_coord_df, zoom=True, zoom_factor=0.7) #it can be helpful to zoom in when checking lomag images
In [72]:
#confirming same is true for high magnification images
input_dict = {'top_left':439,'top_right':440,'bottom_left':512,'bottom_right':511}
TestTileGrid(input_dict, himag_tiles_dir, himag_coord_df, zoom=False, zoom_factor=0)
In [73]:
#update tile grid to reflect most seamless orientation
arranged_grid_lomag = SetTileConfiguration('flipud', untested_grid_lomag)
arranged_grid_himag = SetTileConfiguration('flipud', untested_grid_himag)
In [74]:
#Use updated tile arrangement grids, pre-calculated fiducial points at each magnification, and image dimensions to calculate mapping parameters
P_lomag = Local_to_Global(points_lo, arranged_grid_lomag, [h_lomag, w_lomag])
P_himag = Local_to_Global(points_hi, arranged_grid_himag, [h_himag, w_himag])
mapping_parameters = Fit_By_Points(P_lomag, P_himag, verbose=True)
Verify mapping parameters by mapping a subset of cells and confirming their location at the second magnification¶
In [75]:
df_all_reads = df_lomag
In [76]:
subset_cells = df_all_reads[(df_all_reads["design"]=="ENST00000309318.8_1_3__11") & (df_all_reads["well"]==well_ID)& (df_all_reads["cell_barcode_count_1"]<1)].sample(10)
lomag_subset = np.array(subset_cells[['tile','i_nucleus','j_nucleus']]).astype(int)
In [77]:
#calculating global coordinates from lomag image
calculated_P_lomag_subset = Local_to_Global(lomag_subset, arranged_grid_lomag, [h_lomag, w_lomag])
#using mapping parameters calculated above to calculate global higher mag coordinates
calculated_P_himag_subset = model_TRS(calculated_P_lomag_subset, mapping_parameters, angle='degree')
#use global higher mag coordinates to translate to local higher mag coordinates
reverse_calc_p_himag_subset = Global_to_Local(calculated_P_himag_subset, arranged_grid_himag, [h_himag, w_himag])
In [78]:
#use local lomag coordinates and calculated local higher mag coordinates to look up the same cell at each magnification
#inspect visually that red dot lands on same cell in each image
Plot_Point_Mapping_tif(lomag_subset, reverse_calc_p_himag_subset, lomag_coord_df, himag_coord_df, lomag_tiles_dir, himag_tiles_dir, _vmax=4000)
Use validated mapping parameters to join lower magnification (i.e., genotype) information with higher magnification (i.e., phenotype) information¶
In [79]:
#retrieve local coordinates of cells with reads
df_genotype = df_all_reads[(df_all_reads['cell_barcode_count_0'] > 0) & (df_all_reads["well"]==well_ID)].reset_index(drop=True)
lomag_local_coords = np.array(df_genotype[['tile','i_nucleus','j_nucleus']]).astype(int)
#calculate global coordinate of each cell
calculated_P_lomag = Local_to_Global(lomag_local_coords, arranged_grid_lomag, [h_lomag, w_lomag])
#use mapping parameters validated above to calculate global higher mag coordinates of each cell
calculated_P_himag = model_TRS(calculated_P_lomag, mapping_parameters, angle='degree')
#use global higher mag coordinates to translate to local higher mag coordinates, which are associated with phenotype information
reverse_calc_p_himag = Global_to_Local(calculated_P_himag, arranged_grid_himag, [h_himag, w_himag])
In [80]:
#concatenate genotype and phenotype information
df_highmag_coords = pd.DataFrame(data=reverse_calc_p_himag, columns=['tile_himag','i_himag','j_himag'])
df_genotype[['tile_himag','i_himag','j_himag']] = df_highmag_coords[['tile_himag','i_himag','j_himag']]
df_genotype.to_csv("Well"+str(well_ID)+'_reads_himag-mapped.csv')
In [81]:
df_genotype.sample(10)
Out[81]:
In [ ]:
#the higher mag coordinates can now be used to print out images of cells all unified by some feature described in df_genotype, e.g., sgRNA sequence